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ABSTRACT 

I present an exact and explicit solution to the scalar (Stokes flux intensity) radio 
interferometer imaging equation on a spherical surface which is valid also for non- 
coplanar interferometer configurations. This imaging equation is comparable to vu-term 
imaging algorithms, but by using a spherical rather than a Cartesian formulation this 
term has no special significance. The solution presented also allows direct identification 
of the scalar (spin 0 weighted) spherical harmonics on the sky. The method should 
be of interest for future multi-spacecraft interferometers, wide-held imaging with non- 
coplanar arrays, and CMB spherical harmonic measurements using interferometers. 


1 INTRODUCTION 

Basic radio interferometry deals with narrow fields-of-view 
measured by antenna elements constrained to a plane. Under 
such conditions, i.e. planar brightness distribution and pla- 
nar visibility domain, th e van Cittert-Zernike (vCZ) theorem 
l|Thompson et all 120011 1 states that the brightness and vis¬ 
ibility distributions are two-dimensional Cartesian Fourier 
transforms of each other. An extension of the van Cittert- 
Zernike to a rbitrarily wide fields and non-coplanar arrays 
was given in ICarozzi and Woanl d2009h where it was found 
that the simple Fourier transform relation no longer holds. 

T he genera lized vCZ relation given in 
ICarozzi and Woanl (120091 ) is still similar to the origi¬ 
nal planar vCZ in that the brightness and visibility domains 
are ultimately expressed in Cartesian coordinates. A 
different approach to an interfero metric relation for t he ful l 
celestial sphere was given in lMacphie and Okongw 3 (|l975l) 
which used spherical harmonics in the visibility domain. 
However the main result in that paper was a formula for 
point sources, i.e. the brightness distribution was given 
in terms of delta fu n ctions on the sphere. More recently, 
iMcEwen and Scaifel d2008h used a spherical harmonic 
decomposition of visibility data to obtain the celestial 
sky multipole moments, but their treatment of the radial 
component of the visibility data was not made explicit. 

In what follows I will provide a simple relation, analo¬ 
gous to the vCZ, between a brightness distribution on the 
celestial sphere and its visibility distribution in an arbitrary 
domain — possibly non-coplanar and not necessarily spher¬ 
ical — using a special case of the spherical Fourier-Bessel 
transform rather than using a Cartesian Fourier based trans¬ 
form. 

The vCZ on a sphere relation presented here has sev¬ 
eral practical applications. Spherical harmonics of the sky 
temperature derived from interferometers is of current in¬ 


terest (iKiml 120071: In 3 l200ll) , and in the future there are 
plans for a multi-spacecraft interferometer mission. Such an 
interferometer would observe the full celestial sphere rather 
than a hemisphere, which limits Earth based interferome¬ 
ters. The results are also of interest to observations with 
non-coplanar arrays that currently must deal with the so- 
called ui-term flCornwell and Perlevill992 ), which is a conse¬ 
quence of adapting the two-dimensional, Cartesian Fourier 
transform to work with three-dimensional visibility data to 
produce images of the celestial sphere. Although the imag¬ 
ing technique presented here is naturally suited to extended 
sources and multipole moments, also narrow field-of-view 
interferometers could benefit since, for high dynamic range 
imaging, the trend is to image the entire hemisphere any¬ 
ways in order to handle leakage from beam side-lobes. 


2 A RELATION BETWEEN SKY 

BRIGHTNESS ON CELESTIAL SPHERE 
AND NON-COPLANAR VISIBILITIES 

I start with the scalar intensity component of the extended 
vCZ theorem, i.e. a relation between vi s ibility V and bright¬ 
ness B, as given in ICarozzi and Woanl J2009I ). valid on the 
celestial sphere, which can be written as 

V/(r, k) = j R/(fl fc ) exp (-ik • r) dflfc (1) 

where r is the position vector in the visibility domain, k 
is the wavevector and f= (9k, (pk) are the angular com¬ 
ponents of k on the sphere. The subscripts are used to 
denote that the angles refer to the spherical components of 
the wavevector. Since I will only be concerned with mea¬ 
surements in vacuum, the wavenumber k = |k| is equal the 
frequency used for the visibility measurements divided by 
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the speed of light, uj/c. Note that in eq. 0 the phase refer¬ 
ence position is the origiiQ. The subscript •/ in the equation 
above denotes the Stokes I component, i.e. the scalar flux 
density. In what follows I discard the Stokes I subscript as 
I will only deal with this component. 

The expression 0 actually implies that V fulfills the 
Helmholtz equation, also known as the wave equation. This 
fact is not well appreciated in the radio interferometry lit¬ 
erature, so I present it here. Operating with the Laplace 
operator on eq. 0, one finds that 

V r 2 V + k 2 V = 0 (2) 

which is the Helmholtz equation, or wave equation, in 
the visibility domain. The Helmholtz equation has, besides 
Cartesian solutions, also solutions in spherical coordinates, 
and this suggests that there should be a vCZ relation in 
terms of eigenfun ctions of th e spherical wave equation, which 
are equal to (|jacksonlll999l l 

jt(kr)Y ern (8,cj)), for £ = 0,1,2,...; m =-£,..., £ (3) 

where I invoke the boundary condition that the visibil¬ 
ity should be finite at the origin. Wm(fi) is the standard, 
orthonormal spherical harmonic function with £, m corre¬ 
sponding to the polar and azimuthal quantal number^ re¬ 
spectively, and je (hr) is the spherical Bessel function of the 
first kind. I will call these eigenfunctions spherical wave har¬ 
monics. 

To fully convert eq. 0 into the eigenfunctions given in 
eq. ©, I proceed as fo llows. I use the plane wave decompo¬ 
sition formula, see Ijacksonl d 1999h , 


OO £ 

e lk r = 4n Y (^-i) i jt{kr)Y t . rn (e r ,(j) r )Yl m (ek,(l)k) (4) 

£=0 m=—£ 


where the subscripts » r denote the spherical coordinates of 
the visibility position vector r and r = |r|. When this is 
inserted into eq. 0 it gives 


/ / OO £ \ 

B(n k ) Utt^] Y, (-i) e Mk r )Y ern (0r,<f> r )Y; m {n k )\dn k 

\ £=0 m=—£ ) 

oo £ p 

47T EE (-i ) l jt{kr)Y lm {6r^r) / B(n k )Y e * m (n k )dQ k . 

£=0 m=—£ ^ 

(5) 


Then I expand the brightness distribution into spherical har¬ 
monics 


oo £ 

B(fl k ) = EE (6) 

£=0 m=—£ 

where bim are the multipole moments of the sky. Inserting 


1 That is I have removed the phase reference position so the in¬ 
terferometer is not phased up towards any particular direction. 

2 Since the spherical harmonic quantal numbers £ and m could 
be confused with the standard notation for the direction cosines 
in Cartesian Fourier imaging, the latter will not be used in this 
paper. 


this back into eq. © one obtains 


V = 4.E E (—i) e Jt(kr)Yem(O r ,<l>r) 

£=0 m=—£ 

oo £ 


xf EE 

^ \£=0 m=—£ / 

oo £ 

= 4.E E {-i) e je( kr ) Y em.(9 r ,(t>r)bem (7) 


where I have used the orthogonality relation for the spherical 
harmonic functions 

47T 

U m (n)^* m ,(n)df2 = 8 U '5 mm '- (8) 

Finally I expand the visibility distribution into the spherical 
wave harmonics, eq. ©, with coefficients vim 


oo £ 

v = EE*< mje{kr)Ye m (Q r ). (9) 

£=0 m=—£ 

see Ijacksonl (1 10991 ). Inserting this into the left-hand side of 
eq. 0, one obtains 


oo £ 

^ ^ ^ ^ (^r) 

£=Q m=—£ 

oo £ 

= 4.E E <i(kr)Y em (e r ,(j> r ). (10) 

£=0 m=—£ 

From this equation, due to the orthonormality of the Yg m 
harmonics, one can identify that for any (£, m), 

Vim — 47r( i) bgm. (11) 

Eq. (1111) is an important result, and shows that there 
is a simple proportionality relation between the brightness 
distribution, in terms of 6^ m , and the visibility distribution, 
in terms of vim, with no integration or sum over any domain. 
The simplicity of this result is due to the fact that the spher¬ 
ical harmonic components are eigenfunctions of the mea¬ 
surement equation on the sphere 0 and that these compo¬ 
nents automatically fulfill the Helmholtz dispersion relation 
k 2 = w 2 /c 2 . By contrast, the Cartesian Fourier transform 
consists of plane wave solutions, i.e. point sources, which 
are not eigenfunctions of the measurement equation on the 
sphere and do not automatically fulfill the dispersion rela¬ 
tion which leads to the additional complexity of dealing with 
the w-term, i.e. the third and final wavevector component 
in t he plane wave solutions. 

[McEwen and Scaifel l|2008l ) derived an essentially simi¬ 
lar relationship to (HU. albeit not explicitly. However, they 
did not provide an explicit scheme to derive the harmonic 
coefficients for an arbitrary array. In fact they speculated 
that a stable scheme could be developed, arguing that the 
presence of zeros of the spherical Bessel functions with large 
l would complicate the recovery of the coefficients. I argue 
that the zeros of the spherical Bessel function for some l 
simply mean that that particular £ does not contribute to 
the harmonic coefficient at that point, but e.g. the spherical 
Bessel functions with £ ± 1 will. In the next section I will 
show that the radial part of the visibility can indeed be in¬ 
corporated into the recovery of the spherical harmonics of 
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the sky, and later I will show that, at least for i < 96, it is 
possible to produce images comparable to those made with 
the Cartesian Fourier transform. 


hand side and then get 


poo p'K p2n 00 ^ 

/// EE 

Jo Jo Jo e , =0rn , = _ e , 


vt'm'jf {kr)Y Vm , (9, <t>)je{k 0 r)Yf m (6, tj>) 


3 COMPUTING THE SPHERICAL WAVE 
HARMONIC COEFFICIENTS OF THE 
VISIBILITY DISTRIBUTION 

The result expressed in eq. m shows that spherical har¬ 
monic components of the celestial sky at a some frequency 
w are proportional to a spherical Fourier-Bessel decomposi¬ 
tion of the visibility distribution with the corresponding k. 
Although this vCZ relation is superficially simpler than the 
Cartesian Fourier transform, it still implies a comparable 
computational complexity since the vim components need 
to be determined from the interferometric measurements. 

In practice, an interferometer consists of an array of 
a finite number of antennas from which complex voltages 
are measured, and the visibilities are the complex pow¬ 
ers obtained by cross-correlating between all antenna pairs. 
Thus V can only be sampled at a finite set of Q measure¬ 
ments at points with spherical coordinates which I denote 
as {rt, 9i , Note that from now on I will dispense with 

the r subscripts for the spherical angles in the visibility do¬ 
main that had been used in the previous section. Although 
there are no formal restrictions on the sampling distribution 
for the estimating the vCZ relation in the preceding section, 
certain distributions will be more advantageous than others. 

A detailed discussion of the numerical implementation 
is outside the scope of the present paper, but a direct (non- 
gridded) naive solution for can be derived as follows. 
Consider a visibility dataset Vi(ko) measured in narrow band 
with center frequency too sampled at arbitrary positions. 
These can be seen as a sum of delta functions in the vis¬ 
ibility domain 

Q 

V(k, r, 9, cj>) = ^2 Vi{ko)S(r - r i)5(k — ko) 

i= 1 


00 1 °? 

x r 2 sin 9 drd#d</> = EE ve m / je(kr)j e (k 0 r)r 2 dr 

£=0 m= — £ q 

- 7T<Kfc — fco) ...... 

= /L E ^— ( 14 ) 

e=o m=-e 0 

where I have used the relation 

[ je{kr)j e (k 0 r)r 2 dr m ^ (15) 

Jo 

which is valid for all £, see lLeistedt et all (l2012l l. Integrating 
both the left- and right-hand sides, i.e. the last results of eq. 
m and eq. (USD, over all k and equating these two results 
I find that 

2k 2 Q 

vim{k 0 ) = —■ -y^Vi{ko)je(kori)Yl m (9i,(j)i). (16) 

7T z ‘ 

1=1 

This my main result in terms of providing an explicit, direct 
quadrature rule for computing the spherical wave coefficients 
from arbitrarily placed visibility samples. Note that there is 
no formal restriction on the radial positions of the samples, 
for instance with respect to the zeros of the spherical Bessel 
functions. 

The transform used above to derive eq. m is a type 
of sph e rical Besse l -Four ier decomposition, see lLeistedt et abl 
l|2012lf : iBaddouil (l2010l l. But a crucial difference is that, in 
the present work, the radial component of the wavevector, 
i.e. the wavenumber k, is already known since radio inter¬ 
ferometric visibility data is almost always given as functions 
of frequency in narrow bands, hence the delta function in 
k. Thus the transform is two-dimensional rather than three- 
dimensional for a given frequency. For this reason it may be 
more appropriate to call this something else instead, so I will 
use the term- spherical wave harmonic transform (SWHT). 


= Y, i^e 5{ ~ r - - &)<*(* - M ( 12 ) 

i=1 

where ko = too/c and the factor in the denominator is the 
normalization factor for the delta functions in spherical co¬ 
ordinates. The delta function in the k domain is a simplify¬ 
ing approximation of the spectral density of the frequency 
band response function. Multiplying the right-hand side of 
eq. m with ji(kor)Yl m {9, cj>), and then integrating this over 
a- spherical volume that bounds- the visibility domain re¬ 
sults in 

rrr e 

Jo Jo Jo fc r smd 


4 ALL-SKY IMAGING EXAMPLES 

In this section I apply the SWHT to real radio interferometer 
data to illustrate the imaging technique. I used data from the 
Swedish LOFAR. station, known as SE607, in particular the 
Low Band Array (LBA), which is a ~60 m diameter array 
of 96 crossed dipoles placed in a pseudo-random, co-planar 
pattern covering the frequency range 10 — 90 MHz. 

The dataset I used was a snap-shot, i.e. the cross¬ 
correlations (and auto-correlations) are integrated over a 
short, 10 s, time interval, so that the array can be taken 
to be co-planar. This was chosen since it then can be com¬ 
pared with the ordinary (non-gridded) Cartesian Fourier 
transform. The data was for center frequency 37.1 MHz in a 
192 kHz wide band. I applied eq. GHD to the visibility data 

, \ 2 • Q 1 i/ii i v'G, n \ ■ n -iw* to j. 1 \ an d computed the coefficients up to t = L max = 96, 

x5(k-ko)r sm9drd9d<l>~/SVi(ko)3t(kori)Y em (6i,4>i)5(k-ko). . , , , , ™ ffi . , 

which matches the number of elements, these coefficients 

i= 1 

^g^ere then converted to the sky harmonic coefficients us¬ 
ing eq. (nu, and then these coefficients were used to generate 
For the left-hand side of eq. (USD, 1 insert eq. m and do an image through eq. ©. 

exactly the same the steps as were performed on the right- The result of this SWHT technique is shown in Figure 


© 0000 R.AS, MNRAS 000, 000-000 


















4 T. D. Carozzi 


|T] subplot a). For comparison, subplot b) shows the ordi¬ 
nary Cartesian Fourier transformed image, also known as a 
dirty image, and subplot c) shows a reference model at the 
slightly higher frequency of 50 MHz and with better resolu¬ 
tion. It is clear from this Figure that the SWHT is very sim¬ 
ilar to the Cartesian Fourier transform. The main difference 
is the presence of emissions beyond the telescope horizon, 
i.e. directions apparently below 0 elevation, for instance in 
the South-East corner of subplot b). These emissions are an 
erroneous artifact of the two-dimensional Cartesian Fourier 
imaging technique, since the two Cartesian Fourier compo¬ 
nents, or direction cosines, go from -1 to +1, and there is 
nothing to stop components with absolute value greater than 
1 from contributing to the image. In other words, this illus¬ 
trates the fact that the Cartesian Fourier transform does not 
automatically fulfill the dispersion relation fc 2 = tu 2 /e 2 , as 
already mentioned. 

The run times were slower for the SWHT, but these 
could be improved up by us ing fast spherical h a rmoni c 
transform algorithms, such as iRokhlin and Tygertl ll2006lf . 
which have computation time complexity of the order 
0(N 2 log N), where N is the number of sample points. It 
is possible that an SWHT algorithm could be constructed 
to have a time complexity not much greater than this, mak¬ 
ing it comparable to the w-term imaging algorithms. 


5 CONCLUSIONS 

I have derived a vCZ relation, eq. dm, between a spheri¬ 
cal brightness distribution and an unconstrained visibility 
distribution. I have also presented the SWH transform of 
the visibility data to compute the spherical harmonics of 
the sky from which images can be made. This technique 
was shown to be capable of producing images comparable 
to ordinary dirty images. It should be useful for radio in- 
ferometric imaging of extended sources or for determining 
multipole moments of the celestial sky. It extends naturally 
to visibility data from non-coplanar arrays, and thus the 
technique is comparable to w-term imaging methods. 
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Figure 1. Orthographically projected images of Stokes I flux in arbitrary units of the celestial hemisphere over LOFAR SE607 LB A 
on 3 December 2014, 00:16 UT. The subplot a) image was computed using the spherical Bessel harmonics with max(l) = 96 at 37.1 
MHz. The subplot b) image was computed using a non-gridded two-dimensional Fourier transform also at 37.1 MHz. Subplot c) is an 
reference model image at 50 MHz. The extended emission in the West is the Milky Way. The strong point source in the North-West is 
Cassiopeia A. Emission from the galactic North spur can be seen in the North-East. The circle suggested in these images is the telescope 
local horizon at elevation 0. Note that neither of the images a) or b) have been compensated for the antenna gain pattern. 
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